/*
 * OREKIT-X
 * Copyright 2002-2008 CS Communication & Systemes
 * 
 * Parts of this software package have been licensed to CS
 * Communication & Systemes (CS) under one or more contributor license
 * agreements.  See the NOTICE file distributed with this work for
 * additional information.
 *  
 * This is an experimental copy of OREKIT from www.orekit.org.
 * Please use the original OREKIT from orekit.org for normal work
 * unrelated to this research project.
 * 
 * Licensed under the Apache License, Version 2.0 (the "License"); you
 * may not use this file except in compliance with the License.  You
 * may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
 * implied.  See the License for the specific language governing
 * permissions and limitations under the License.
 */
package ore;

import java.io.Serializable;

import org.apache.commons.math.geometry.Vector3D;
import ore.time.AbsoluteDate;
import ore.time.TimeStamped;
import ore.PositionVelocity;


/**
 * Position-Velocity model based on Chebyshev polynomials.  This class
 * represents the elements of the piecewise ephemerides for solar
 * system bodies like JPL DE 405 ephemerides.
 * 
 * @author Luc Maisonobe
 */
public class PositionVelocityChebyshev
    extends Object
    implements TimeStamped, Serializable
{

    public final AbsoluteDate start;
    public final double duration;

    private final double[] xCoeffs;
    private final double[] yCoeffs;
    private final double[] zCoeffs;


    public PositionVelocityChebyshev(AbsoluteDate start, double duration,
                                     double[] xCoeffs, double[] yCoeffs, double[] zCoeffs)
    {
        super();
        this.start    = start;
        this.duration = duration;
        this.xCoeffs  = xCoeffs;
        this.yCoeffs  = yCoeffs;
        this.zCoeffs  = zCoeffs;
    }


    /** {@inheritDoc}
     */
    public final AbsoluteDate getDate() {
        return this.start;
    }
    /**
     * @return Model validity duration in seconds
     */
    public final double getValidityDuration() {
        return this.duration;
    }
    /**
     * @param date Check
     * @return Date is in validity range
     */
    public final boolean inRange(final AbsoluteDate date) {
        final double dt = date.durationFrom(start);
        return (dt >= 0) && (dt <= this.duration);
    }
    /**
     * @param date Date for position-velocity 
     * @return Position-velocity at specified date
     */
    public final PositionVelocity getPositionVelocity(final AbsoluteDate date) {
        /*
         * Normalize date
         */
        final double t = (2 * date.durationFrom(start) - duration) / duration;
        final double twoT = 2 * t;
        /*
         * Initialize position 
         */
        double pKm1 = 1.0;
        double pK   = t;
        double xP   = xCoeffs[0];
        double yP   = yCoeffs[0];
        double zP   = zCoeffs[0];
        /*
         * Initialize derivatives
         */
        double qKm1 = 0.0;
        double qK   = 1.0;
        double xV   = 0.0;
        double yV   = 0.0;
        double zV   = 0.0;
        /*
         * Combine polynomials by applying coefficients
         */
        for (int k = 1, z = xCoeffs.length; k < z; ++k) {
            /*
             * Position
             */
            xP += xCoeffs[k] * pK;
            yP += yCoeffs[k] * pK;
            zP += zCoeffs[k] * pK;
            /*
             * Velocity
             */
            xV += xCoeffs[k] * qK;
            yV += yCoeffs[k] * qK;
            zV += zCoeffs[k] * qK;
            /*
             * Next polynomial 
             */
            {
                final double pKm2 = pKm1;
                pKm1 = pK;
                pK   = twoT * pKm1 - pKm2;
            }
            /*
             * Next derivative
             */
            {
                final double qKm2 = qKm1;
                qKm1 = qK;
                qK   = twoT * qKm1 + 2.0 * pKm1 - qKm2;
            }
        }

        final double vScale = (2.0 / duration);
        return new PositionVelocity(new Vector3D(xP, yP, zP),
                                    new Vector3D(xV * vScale, yV * vScale, zV * vScale));
    }
    public final int compareTo(TimeStamped stamped) {
        return this.getDate().compareTo(stamped);
    }
}
